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A  NEW  METHODOLOGY  FOR  THE  NUMERICAL 
SIMULATION  OF  STRAIN  SOFTENING 
IN  INELASTIC  SOLIDS 


J.  C.  SIMO 

Division  of  Applied  Mechanics 
Department  of  Mechanical  Engineering 
Stanford  University,  Stanford,  CA  94305 


Abstract 

This  document  summarizes  the  work  performed  by  the  author  and  his  associates  at  Stanford 
on  a  new  approach  to  the  analysis  and  simulation  of  localization  arising  in  inelastic  solids  that 
exhibit  strain-softening  response.  In  addition,  the  document  describes  new  results  pertaining  to 
the  extension  of  these  ideas  to  the  finite  deformation  regime.  Such  an  extension  will  be  exploited 
numerically  in  follow-up  publications.  The  techniques  described  below  lead  to  the  systematic  con¬ 
struction  of  numerical  methods  that  completely  eliminate  the  strong  mesh  dependence  exhibited 
by  conventional  treatments  of  the  problem. 


1.  Introduction. 

In  a  recent  paper  Simo,  Oliver  &  Armero  [1993]  identified  key  qualitative  features 
exhibited  by  the  response  of  rate-independent  inelastic  solids  in  the  presence  of  strain 
softening  and  demonstrated  that  these  features  are  consistent  with  solutions  possessing  a 
discontinuous  displacement  field  in  the  (quasi-static)  rate-independent  theory.  The  analy¬ 
sis  performed  in  this  reference  for  the  general  three-dimensional  problem  shows  that  strong 
discontinuities  (i.e.,  jumps  in  the  displacement  field)  are  consistent  with  rate-independent 
softening  response  in  the  isothermal  quasi-static  regime,  provided  that  the  softening  modu¬ 
lus  is  re-interpreted  in  a  distributional  sense.  In  particular,  for  linear  softening  the  inverse 
of  the  plastic  softening  modulus  becomes  a  constant  times  a  delta-function  localized  at  the 
discontinuity  if  the  problem  is  to  remain  mathematically  well-posed.  In  a  subsequent  pa¬ 
per,  Oliver  &  Simo  [1994]  extended  the  preceding  analysis  to  a  widely  used  class  of  isotropic 
damage  models  and  showed  that  the  same  conclusions  remain  valid.  In  particular,  we  have: 

i.  The  distributional  re-interpretation  of  the  softening  law  leads  to  a  compelling  re¬ 
interpretation  of  the  softening  modulus,  related  to  the  fracture  energy-released  rate 
in  the  solid. 
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ii.  Displacement  solutions  in  the  presence  of  strain-softening  necessarily  yield  discontin¬ 
uous  displacement  fields  and,  therefore,  strain  fields  which  exhibit  delta-measures. 

iii.  Any  numerical  methods  suitable  for  strain  localization  must  be  able  to  reproduce 
strain  fields  with  delta-measure.  This  requirement  together  with  condition  i.  yield 
numerical  solutions  that  exhibit  no  mesh- dependence  whatsoever. 

Strong  discontinuities  are  understood  in  this  work  in  the  sense  of  solutions  exhibiting 
jumps  in  the  displacement/velocity  field  across  material  surfaces  (or  sets  of  zero  measure). 
The  displacement  gradients  on  both  sides  of  the  surface  of  discontinuity  then  differ  by  a 
bounded  measure  (a  delta-function)  while  the  traction  vector  remains  continuous.  The 
classical  example  afforded  by  the  slip-lines  of  rigid  plasticity,  i.e.,  bands  of  zero  thickness 
and  localized  plastic  deformation.  By  contrast,  solutions  exhibiting  weak  discontinuities 
remain  continuous,  while  the  displacement  gradient  (not  the  displacement  itself)  experi¬ 
ences  a  jump  across  the  surface  of  discontinuity.  A  typical  example  is  provided  by  the 
weak  shocks  observed  in  elastic  materials  exhibiting  phase  transitions.  In  spite  of  some 
confusion  in  the  literature,  the  numerical  analysis  issues  involved  in  the  simulation  of  weak 
shock  are  significantly  different  from  those  arising  in  the  simulation  of  strong  shocks.  The 
former  case  is  treated  in  reference  Mamiya  k  Simo  [1994]. 

A  number  of  numerical  techniques  have  been  proposed  in  an  effort  to  achieve  high 
resolution  of  strong  discontinuities,  the  typical  example  being  the  discontinuous  Galerkin 
finite  element,  method.  Unfortunately,  these  techniques  are  often  cumbersome  to  imple¬ 
ment,  rather  expensive,  and  do  not  fit  the  framework  for  developing  high  resolution  schemes 
for  sharp  shock  capturing.  In  Simo,  Oliver  k  Armero  [1993]  and  Simo  k  Oliver  [1994], 
we  have  shown  that  the  Assumed  Enhanced  Strain  (AES)  method  of  Simo  k  Rifai  [1990] 
provides  an  ideal  framework  for  developing  high  resolution  schemes  for  sharp  shock  cap¬ 
turing.  The  AES  methodology  has  been  analyzed  from  a  mathematical  point  of  view  in 
Reddy  k  Simo  [1992]  and  extended  to  the  finite  deformation  regime  in  Simo  k  Armero 
[1992]  and  Simo,  Armero  k  Taylor  [1993].  The  underlying  idea  for  the  problem  at  hand  is 
to  replace  the  strain  field  by  an  enhanced  strain  field  designed  to  replicate  delta-functions 
within  a  typical  element.  The  explicit  construction  is  described  in  detail  in  Simo  k  Oliver 
[1994]  within  the  context  of  the  infinitesimal  theory  and  below  for  the  full  finite  deforma¬ 
tion  problem.  A  similar  idea  also  works  for  the  motion  of  weak  shocks  in  phase  transition 
problems  (see  Mamiya  k  Simo  [1994]). 

An  outline  of  the  rest  of  this  report  is  as  follows.  Section  2  generalizes  the  basic  results 
described  in  Simo  k  Oliver  [1994]  on  the  kinematics  and  balance  laws  for  the  problems 
with  strong  discontinuities  to  the  full  finite  deformation  regime.  Section  3  presents  com¬ 
pletely  new  results  on  finite  strain  plasticity  in  the  presence  of  strong  discontinuities  that 
generalizes  to  the  finite  strain  regime  the  techniques  described  in  Simo,  Oliver  k  Armero 
[1993]  and  summarized  above.  Finally,  Section  4  describes  a  new  finite  element  method  for 
the  accurate  resolution  of  strong  discontinuities  that  generalizes  to  the  finite  strain  regime 
the  technique  described  in  Simo  k  Oliver  [1994].  In  sharp  contrast  with  recent  approaches, 


J.C.  Simo 


3 


see  e.g.  Larsson,  Runesson  &  Ottosen  [1994],  in  the  techniques  developed  in  the  course  of 
this  research  the  shock  does  not  need  to  coincide  with  element  sides. 

As  part  of  the  proposed  research  effort,  we  have  developed  new  constitutive  models 
for  geomaterials  in  Simo  &  Meschke  [1993].  In  addition,  a  rather  attractive  material  model 
for  plane  concrete  that  incorporates  anisotropic  damage  is  presented  in  Govindjee,  Kay  & 
Simo  [1994].  Lack  of  space  prevents  us  to  describe  in  detail  these  developments.  Similarly, 
no  details  will  be  given  on  either  the  mathematical  analysis  of  enhanced  strain  meth¬ 
ods  presented  in  Reddy  &  Simo  [1992]  or  the  class  of  assumed  enhanced  strain  methods 
described  in  Simo,  Armero  &  Taylor  [1993]. 


2.  Kinematics  and  Balance  Laws  in  the  Presence  of  shocks. 

We  summarize  below  the  basic  notation  employed  throughout  this  work  and  introduce 
a  fundamental  decomposition  of  the  displacement  field  into  a  regular  part  and  a  discontin¬ 
uous  part.  Such  a  decomposition  plays  a  key  role  in  proposed  numerical  treatment  of  the 
problem.  We  conclude  this  section  with  the  statement  of  the  weak  form  of  the  equations 
and  a  rigorous  derivation  of  the  local  equilibrium  equations. 


FIGURE  2.1.  Material  surface  5  in  the  domain  12,  with  unit  normal  N  : 
S  —*  S2,  and  normal  paramaterization  X  =  Y  +  tjN  with  i?  €  [—  h,h].  The 
surface  5  divides  both  f2  and  C2h  C  12  into  two  subdomains  labeled  f2±  and 
S7^_  C  J2±,  respectively. 


2.1.  Spatially  Discontinuous  Motions.  Basic  Notation. 

We  denote  by  i?  C  Rndim  the  reference  configuration  of  a  continuum  body  with  smooth 
boundary  dQ  and  particles  labeled  by  X  €  Q.  The  starting  point  of  the  analysis  of 
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solutions  exhibiting  strain  softening  is  a  motion  p  :  Q  x  I  — *  R"dim,  which  is  assumed  to 
experienced  a  jump  discontinuity  \p\  ^Oona  smooth  material  surface  S  C  Q.  Here  ! 
denotes  the  time  interval  of  interest.  Since  S  is  a  material  surface  we  have  S  =  0.  Points 
in  S  are  denoted  by  Y,  with  the  notation  indicated  in  Figure  1  assumed  throughout.  For 
simplicity  we  assume  that  S  is  smooth,  with  unit  normal  vector  field  denoted  by  N(Y )  so 
that 

S=[y  =  Y(Z\t2)  :  (£\£2)  €  B }  ,  (2.1a) 

where  Y  :  B  — *  Q  is  a  smooth  global  parametrization.  The  unit  normal  field  to  S  is  then 
given  by 

jv  =  $(?,?)  —  v,i  xy,2  /  nf.j  xY,2  ||.  (2.1b) 

A  parametrization  of  S  induces  in  turn  a  convenient  parameterization  valid  in  a  neigh¬ 
borhood  Qh  =  S  x  (—h,  h )  of  S,  known  as  a  normal  parametrization  and  defined  via  the 
formula 

X(^^2,n)  +  for  —  h  <  rj  <  h.  (2.2) 

The  boundary  dQh  of  the  set  Qh  is  composed  of  the  surfaces  shown  in  2.1,  so  that 

8Qh  =Sh_VSl  Ui^Uit,  (2.3) 

where  rh  =  Tt  U  is  the  intersections  of  Qh  with  the  outer  boundary  dQ  of  Q.  Given 
an  arbitrary  function  ip  :  Qh  x  I  — >  R,  we  let 

(2-4«) 


With  a  slight  abuse  in  notation,  we  shall  use  the  same  symbol  to  denote  both  the  function 
ip  and  its  representative  ip  defined  by  (2.4a)  in  the  normal  parametrization.  In  addition,  we 
denote  by  {G1,  G2,  iV}  the  basis  dual  to  the  curvilinear  basis  {X,i  ,  X,2  ,N}  associated 
with  the  normal  parametrization,  so  that  Xa  •  —  8 @  for  a,  ^  =  1,  2.  It  follows  that  the 

gradient  of  ip  can  be  written  as 

GRAD  [ip]  =  GRADs  [ip]  +  [dip  /  dr)}  N ,  (2.46) 


where 


*Mim  1 

GRADSW~  Y,  (W/d(°}G°, 


(2.4c) 


CK=  1 

Clearly,  by  construction  we  have  N  ■  GRADs [t/>]  =  0.  Finally,  the  motion  tp(X,t)  is 
subjected  to  the  usual  essential  and  boundary  conditions.  Explicitly,  we  assume  that  on 
parts  ru  C  dQ  and  P*  Cf  we  have 


ip  —  g  on  ru  x  I  and  Pv  =  t  on  P*  x  5, 


(2.5) 


where  g  and  t  are  prescribed  vector  fields  on  <9 12,  u  is  the  unit  outward  normal  field  to  d£2 
and  P  denotes  the  nominal  stress  tensor.  We  assume  that  P„  UPt  =  dQ  while  Pu  D  P*  =  0. 
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2.2  Decomposition  of  the  Spatially  Discontinuous  Motion 

The  numerical  approach  described  below  relies  crucially  on  a  decomposition  of  the 
motion  into  a  regular  or  smooth  part  and  a  discontinuous  part.  The  regular  part  in  such 
a  decomposition  is  assumed  to  obey  the  same  essential  boundary  conditions  (2.5)i  as  the 
motion.  Accordingly,  we  set 

<p(X,t)=  +Ms(X)M(Y,tl  for  (X,t)  £0x1.  (2.6) 

Regular  Part  Discontinous  Part 

Here,  [^?](y,t)  is  the  jump  experience  by  (p  across  S  defined  in  the  usual  fashion,  see 
e.g.,  Truesdell  &  Toupin  [1960],  and  Ms  :  &h  -»  R  is  a  prescribed  function  subject  to  the 
following  two  conditions: 

i.  The  jump  across  S  is  [Ms]  =  1. 

ii.  The  support  of  Ms  is  suppfMs]  =  fih  with  Ms  =  0  on  S±. 

A  convenient  expression  for  the  function  Ms  is  obtained  by  means  of  the  following  con¬ 
struction.  Relative  to  a  normal  parametrization  X  —  Y  +  rj  N  we  set 

MS(X)  =  Hs{ri)  -  i>h(X)  where  Hs(v)  =  j  J  -f  ^  <  o’  (2J) 

is  the  Heaviside  function.  The  function  ijjh  is  smooth  and  completely  arbitrary,  except  for 
the  following  two  requirements: 

<Ph(X)  =  0  for  X  €  Q-\nh_  and  $\X)  =  1  for  X  €  (2.8) 

The  properties  of  the  Heaviside  function  imply  that  condition  i  holds  for  Ms  defined  by 
(2.7),  while  restriction  (2.8)  ensures  that  condition  ii  also  holds. 

Remark  2.1.  For  the  case  n<nm  =  15  collapses  to  a  point,  say  x  =  £,  and  definition  (2.7) 
reduces  to  Ms(rj)  =  +v)~  +  7?)-  Tiie  restriction  (2.8)  on  iph  merely  becomes 

+  rj)  =  0  for  r)  <—h  and  +  rf)  =  1  for  r]  >  h.  (2.9) 

■  The  derivative  of  Ms  is  obviously  given  by  M's{rj)  =  +  r?)  -  ?phl(£  +  rj ).  ■ 

For  subsequent  use,  we  recall  that  GRAD  [if  s]  =  8s  N  for  ndim  >  1?  where  8s  is 
the  Dirac  delta  measure  on  the  surface  S.  Using  this  elementary  result  in  the  theory  of 
distributions,  the  derivative  of  the  function  Ms  takes  the  form 

GRAD  [Ms]  =  8SN  -  GRAD[</>fe]. 


(2.10) 
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We  remark  that  the  smoothness  restriction  on  the  regular  part  <p  of  the  motion  is  intro¬ 
duced  only  for  simplicity.  Continuity  of  <p  is  all  that  is  required  for  the  arguments  below 
to  remain  valid. 


2.3  The  Weak  and  Local  Form  of  the  Equilibrium  Equations 

Again  for  simplicity,  attention  will  be  restricted  in  the  developments  below  to  the 
quasi-static  problem  under  the  assumption  of  small  deformations.  Our  goal  is  a  careful 
statement  of  the  local  equations  of  equilibrium  for  the  problem  at  hand.  Our  point  of 
departure  is  the  weak  form  of  the  equations  as  defined  by  the  classical  principle  of  virtual 
work,  i.e., 

f  P  :  GRAD [17]  dQ  =  [  f  rjdQ  +  f  t  rjdr ,  (2.11) 

Jo  Jo  J  rt 

for  all  admissible  test  functions  (virtual  displacements)  rj  6  V,  where  f  stands  for  the 
prescribed  body  force.  Consistent  with  the  form  (2.6)  adopted  for  the  displacement  field, 
the  space  V  of  admissible  variations  is  defined  by 


:=  j*7  =  V 


+  Ms  (3  :  rj 


ru 


=  0  and  /3  :  S  — ►  IRndim  is  arbitrary!  ,  (2.12) 


J 


where  the  regular  part  rj  of  rj  €  V  is  smooth.  We  have  the  following  result. 


Lemma  2.1.  The  weak  form  (2.11)  along  with  (2.12)  yield  the  local  equilibrium  equation 


div[P]  +  /  =  0  in  f2\S  x  I 
Pv  =  t  on  rt  x  I 


(2.13a) 


supplemented  with  the  following  two  conditions: 

| PNJ  =  { P+  -  P-  }  N  =  0  and  Ps  N  =  P+  N.  (2.136) 

Here  Ps  N  denotes  the  traction  vector  acting  on  the  surface  S. 

Proof:  We  proceed  in  two  steps,  (i)  First,  choose  regular  test  functions  with  (3  =  0,  so 
that  GRAD [17]  =  GRAD [fj]  with  rj  arbitrary.  Integration  by  parts  then  yields 

f  P  :  GRADfa]  dQ  =  —  /  fj  DIV[P]  dQ 

Jn  Jo 

+  J  - dS+ Jrrr(PN)dr.  (2.14) 

By  inserting  (2.14)  into  (2.11)  a  standard  argument  yields  (2.13a)  and  (2.13b)i. 
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(ii)  Second,  choose  test  functions  with  the  regular  part  rj  =  0  and  function  Ms  defined 
by  (2.7).  From.  (2.10)  we  have 

GRAD [77]  =  -P®  GRAD[^A]  +  MS  +  6SP  ®  N  +  GRAD 5 [/3],  (2.15) 

where  GRAD5[/3]  =  P,a  see  (2.4).  The  contribution  of  the  first  term  in 

(2.15)  to  the  weak  form  (2.11)  is  evaluated  as  follows.  Integration  by  parts  along  with  the 
local  equilibrium  equation  (2.13a)i  gives 


®  GRAD[V>&])  dQ 


[  m\[xl>hPf5]dQ 

Jnh 

-  [  iph[P-  DIV[P]  +  p  :  GRADs [/3]]  dQ 

f  ij)hp  ■  PN  dS  -  f  iphp-PNdS 
J  «S_|_  J  5- 

-  J  $kP-\P\NdS 

+  [  iph[f  ■  P  -  P  :GRABs[P]\  dQ.  (2.16) 

Jnh 


Use  of  the  boundary  conditions  ij)h  =  1  on  and  —  0  on  S'L,  together  with  the  jump 
condition  (2.13b)2  and  expression  (2.7)  for  Ms,  yields 


GRADED  dQ  =  [  p-PNdS-  I  [ f  P-P :  grad5[/?]]  Ms  dQ 
Js  Jnh 

+  f  [f  p-P  :  GRADs [/3]]  dQ.  (2.17) 

Jn^ 


Inserting  (2.15)  into  the  weak  form  (2.11)  of  the  equilibrium  equations  and  making  use  of 
(2.17)  gives  the  following  result: 

f  P-PsNdS-  f  p-PNdS-  f  [f  ■  P  -  P  :  GRAB S[P]]  dQ  =  0,  (2.18) 

Js  Js £  Jnh+ 

for  any  P  :  S  — >  R”dim.  The  last  step  in  the  proof  involves  the  evaluation  of  the  integral 
involving  the  term  P  :  GRAD5[/3]  =  P,a  Pga.  Since  in  a  normal  paxametrization  Q\  = 
S  x  [0,  h],  using  integration  by  parts  and  the  assumption  that  the  lateral  boundary  = 
ph  y  ph  js  traction  free  it  is  easily  shown  that 


[  P -.GRAB  s[P}dQ  =  -  f  p-Blv[P}dQ  +  f  P -(PN)^  dSdrj.  (2.19) 
Jab  Jn*  Jsx[o,h] 
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Inserting  the  local  equilibrium  equation  (2.13a)i  into  (2.19)  yields 

-  f  [ f  (3- P :  grad5[/3]]  dQ  =  [  (3-  ( PN )  dS  -  [  (3  ■  ( P+N )  dS.  (2.20) 

Jn\  Js 

By  combining  (2.20)  and  (2.18)  we  conclude  that  the  weak  form  of  the  equilibrium  equa¬ 
tions  finally  reduces  to 


J  f 3  ■  [PSN  -  P+N }  dS ,  V/3  :  <S 


P^dim 


(2.21) 


A  standard  argument  then  yields  result  (2. 13b) 2-  B 

Remark  2.2.  Results  (2.13a)  and  (2.13b)x  are  classical  for  problems  exhibiting  discon¬ 
tinuities.  Result  (2.13b)2,  on  the  other  hand,  is  non-conventional  and  arises  from  the 
decomposition  (2.6)  of  the  displacement  field.  The  different  physical  significance  of  these 
two  conditions  becomes  clear  if  we  replace  the  shock  S  by  the  neighborhood  fih  with  finite 
thickness  h.  Condition  (2.13b)x  then  becomes  a  restatement  of  continuity  of  the  traction 
vector  across  the  material  shock  <S,  in  the  sense  that  the  traction  vectors  on  <S_£  and  <S_ 
coincide,  while  condition  (2.13b)2  is  the  assertion  that  traction  vector  on  the  shock  itself 
S  equals  the  traction  vector  on  S+  ahead  the  shock.  B 


3  Finite  Plasticity  in  the  Presence  of  Softening 

Classical  rate  independent  models  of  plasticity  lead  to  ill— posed  initial  boundary  value 
problems  in  the  presence  of  strain— softening  (as  opposed  to  strain— hardening)  response. 
This  lack  of  well-posedness  manifests  itself  in  a  computational  setting  in  a  strong-mesh 
dependence  of  the  numerical  solution.  Within  the  scope  of  the  infinitesimal  theory,  it  was 
shown  in  Simo,  Oliver  &  Armero  [1993]  that  well— posedness  in  strain— softening  models  of 
rate-independent  plasticity  is  restored  if  the  hardening  law  is  reinterpreted  in  a  distribu¬ 
tional  sense.  Moreover,  this  reformulation  of  the  Hardening  law  lends  itself  to  a  compelling 
physical  interpretation.  The  goal  of  this  section  is  to  generalize  these  results  to  the  full 
finite  deformation  theory. 

3.1  Deformation  Gradients  and  Spatial  Rates  of  Deformation 

The  first  step  in  the  analysis  is  to  compute  the  deformation  gradient  and  the  associated 
rates  of  deformation  tensors  for  a  discontinuous  motion  of  the  form  (2.6).  For  convenience, 
Ms  is  assumed  to  be  given  by  (2.7).  Accordingly,  setting 

F  :=  GRAD[^]  +  jWsGRADs  -  {p}  ®  GRAD^], 


(3.1) 
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and  recalling  that  GRAD[iJs]  =  SsN,  the  deformation  gradient  F  :=  GRAD[<£>]  can  be 
written  as 

F=F  +  8sl<p}®N  in  f2xl.  (3.2 a) 

Alternatively,  by  defining  the  material  jump  J  as 

J  =  F~1lcp  J,  '  (3.3) 

expression  (3.2)  for  the  deformation  gradient  can  be  written  in  the  following  entirely  equiv¬ 
alent  form: 

F  =  FF*  where  F*  :=  1  +  6s  J  0  N.  (3.2 b) 

This  kinematic  decomposition,  which  arises  in  the  present  context  without  any  a-priori 
kinematic  assumption,  is  formally  identical  to  the  multiplicative  decomposition  F  =  FeFp 
of  the  deformation  gradient  postulated  at  the  outset  in  certain  formulations  of  plasticity, 
see  e.g.,  Lee  [1969],  Mandel  [1974]  and  Simo  [1992,  1994]  among  others.  For  the  problem 
at  hand,  the  elastic  part  is  JFe  =  F  while  the  plastic  part  becomes  Fp  =  F*.  In  order  to 
circumvent  technical  difficulties,  the  following  additional  hypothesis  is  introduced. 

ASSUMPTION  3.1.  The  unit  vector  N  normal  to  the  material  discontinuity  S  is  orthogonal 
to  the  material  jump  J,  so  that  the  following  equivalent  conditions  are  assumed  to  hold: 

N  -  J  =  0  4=^  n  •  |y>]  =  0  where  n  :=  F~TN.  (3-4) 

It  follows  from  (3.4)2  that  n  can  by  interpreted  as  the  vector  normal  to  the  surface  <p(S) 
in  the  current  configuration  of  the  solid.  ■ 

Expressed  in  physical  terms,  condition  (3.4)  asserts  that  only  slips  (i.e.,  jumps  \<p\ 
which  are  tangential  to  S )  are  allowed.  Under  such  a  restriction,  the  inverse  of  the  defor¬ 
mation  gradient  F*  takes  the  simple  form 

F*"1  =  1  -  6S  J  0  JV;  (3.5) 

a  result  which  is  immediately  verified  by  a  direct  computation. 

Remark  3.1.  If  the  normal  n  to  <p(S)  is  not  orthogonal  to  the  jump  [<£>],  it  is  easily  verified 
that  expression  (3.5)  remains  formally  valid,  provided  that  the  delta  measure  6S  replaced 
by  the  factor  [1  —  •  n]-1^.  H 

The  next  step  is  to  compute  the  spatial  velocity  gradients  l  :=  GRADf^F-1  and 
l  :=  GRAD[^]F_1  associated  with  the  motions  <p  and  (p,  respectively.  A  straightforward 
manipulation  then  yields  the  following  standard  expressions  in  continuum  mechanics: 

Z  =  FF-J  and  T=P'F~1.  (3.6) 


According  to  the  preceding  interpretation,  we  can  view  l  and  l  as  the  total  and  the  elastic 
spatial  velocity  gradients,  respectively.  The  following  result  determines  the  plastic  velocity 
gradient. 

LEMMA  3.1.  The  local  velocity  gradient  l*  associated  with  the  part  F*  of  the  deformation 
gradient  obeys  the  additive  decomposition 

l  =  1  +  1*  with  l*  6s  [<£>!)  ®  (3.7a) 

and 

Lvlvl  :=  ~  where  ®  :=  0  £_1  (3‘76) 

is  the  Lie  derivative  relative  to  the  spatial  velocity  field  v. 

PROOF.  Time  differentiation  of  the  relation  F  =  FF*  yield,  after  post-multiplication  by 
F -1,  the  result 

dF* 

l  =  T+  F  [L*]  F-1  where  L*  :=  —  F*~\  (3.8) 

Now,  time  differentiation  of  expression  (3.2b)2  and  use  of  result  (3.5)  yields 

L*  =  \SSJ  (8»  JV]  [1  -  6s  J  ®N]  =  6S  J  O  N,  (3.9) 

since  AT  •  J  =  0  by  assumption  (3.4).  Therefore,  we  have 

F  [j L *]  F-1  =  6S(FJ)  ®n  =  6s  L9[<p }  ®  n.  (3.10) 

and  the  result  follows  by  inserting  (3.10)  into  (3.8).  ■ 

REMARK  3.2.  A  standard  result  in  continuum  mechanics  shows  that  the  Lie  derivative 
Lel<p]  is  objective.  Since  the  normal  n  is  also  objective,  it  follows  that  the  full  spatial 
velocity  gradient  gradient  l*  is  objective  in  the  standard  sense,  i.e., 

l*+  =  Ql*QT  if  x^x+  :=r  +  Qx,  (3.11) 

for  arbitrary  translations  re  R3  and  arbitrary  rotations  Q  €  50(3).  On  the  other 
hand,  only  the  symmetric  parts  d  :=  sym[Z]  and  d  sym[Fj  — the  so-called  spatial  rate 
of  deformation  tensors —  are  objective.  An  identical  result  holds  for  finite  deformation 
multiplicative  plasticity,  see  Simo  [1994].  We  shall  denote  by  w  :=  skew[Z],  w  =  skew[Z] 
and  w*  =  skew[/*]  the  spin  tensors.  ■ 
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3.2  Analysis  of  Discontinuous  Solutions.  Localization  Conditions 

Consider  next  a  conventional  model  of  classical  rate— independent  plasticity  with  isotropic 
hardening/softening  law  expressed  in  rate  form  via  the  following  constitutive  equations: 

T  =  C  [l  —  P] 

P  =  XdT<f>(r,q) 
a  =  \dq<j>(T,q) 

A  >  0,  4>(t,  q)  <  0  and  A  <f(r,  q)  =  0. 

Here  r  :=  PF~T  is  the  Kirchhoff  stress  tensor  and  r  :=  f  -  wr  +  rw  its  Jaumann 
derivative,  <^>(t,  q)  :=  <^»(t)  ~t~  3  —  &Y  is  the  yield  criterion  and  cry  >  0  is  the  flow  stress.  We 
assume  that  $(■)  is  convex  function,  homogeneous  of  degree  one.  In  addition,  the  spatial 
elasticity  tensor  c  is  assumed  to  possesses  the  usual  symmetry  conditions  and  remains 
positive  definite,  in  the  sense  that  t  •  c-1t  >  A>  1 1  'r  j  | 2  with  flo  >  0  and  arbitrary  t  =  t  . 
To  illustrate  the  key  results,  it  suffices  to  consider  the  simplest  model  of  softening  plasticity 
in  which 

a  =  — 7f-1  q  with  =  constant  <  0.  (3.13) 

The  case  Ti  >  0  corresponds  to  hardening  plasticity.  Observe  that  the  flow  rule  (3.12)2 
implies  that  the  plastic  spin  wv  =  0  since  dr<f>  is  symmetric. 

REMARK  3.3.  The  preceding  model  can  be  shown  to  arise  as  the  rate  form  of  multiplicative 
multiplicative  plasticity,  with 

le  =  jpejpe-i  and  ip  ■-  FeLpFe~\  where  LP  =  FPFP~1.  (3.14) 

The  Kirchhoff  stress  tensor  is  defined  in  terms  of  a  stored  energy  function  W (Ce)  via  the 
hyperelastic  constitutive  equation 

r  =  Fe  [2 dc*W(Ce)\  FeT  where  Ce  =  FeTFe.  (3.15) 

The  spatial  elastic  moduli  c  are  given  by  the  standard  relation 

d2W 

=  F&Fh*kK*TL  <3-16) 

Finally,  the  local  plastic  dissipation  T>  predicted  by  model,  defined  as  he  local  stress  power 
minus  the  rate  of  change  of  the  internal  energy,  equals  the  flow  stress  times  the  plastic 
slip,  i.e., 


(3.12) 


V  =  Ary  >  0,  since  A  >  0. 


(3.17). 
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A  detailed  proof  of  the  preceding  results  is  given  in  Simo  [1994].  ■ 

The  goal  of  the  analysis  below  is  to  identify  the  conditions  that  render  the  constitutive 
model  (3.12)  with  the  softening  law  (3.13)  mathematically  consistent  with  expression  (3.7a) 
for  the  spatial  velocity  gradient  Z,  which  involves  a  delta- measure.  We  proceed  in  two  steps. 

i.  First,  by  inserting  (3.7a, b)  into  into  the  rate  constitutive  equations  (3.12)  one 
arrives  at  the  following  expression 

c_1r  —  f  =  5ssym  <8>  «•]  —  A dT<f> ■  (3.18) 

Now,  the  left-hand-side  of  this  equation  is  a  piecewise  smooth  function  since  T  is  piecewise 
smooth  and  the  stress  field  r  can  experience  at  most  jump  discontinuities  (i.e.,  r  is  of 
bounded  oscillation).  The  right -hand-side  of  this  equation,  on  the  other  hand,  exhibits 
a  term  involving  a  ^-measure.  Since  both  factor  multiplying  8$  is  a  smooth  function 
defined  on  S  and  the  term  dr<f>  are  piecewise  smooth,  equation  (3.18)  can  make  sense 
mathematically  only  if  the  plastic  multiplier  A  itself  is  a  distribution  of  the  form 

A  =  A  +  A 8s  where  A  :  S  — >  R  (3.19) 

is  smooth  and  non-negative  and  corresponds  to  localized  plastic  deformation  on  the  surface 
S.  The  function  A  :  Q  — >  R  is  a  piecewise  smooth  and  non-negative  and  corresponds  to 
diffuse  yielding  in  the  solid.  Relations  (3.18)  and  (3.19)  then  imply  that 

sym  [v?]  ®  n]  =  A dT<j>  on  5x1,  (3.20) 

We  shall  refer  to  this  result  as  the  localization  condition.  Without  loss  of  generality,  in 
what  follows  it  will  be  assumed  that  A  =  0  (localized  plastic  deformation)  so  that  (3.18) 
reduces  to 

t  =  cT  in  {2\S  x  8  and  A  =  X8s-  (3.21) 

In  other  words,  the  region  !2\<S  remains  elastic,  in  agreement  with  the  simplifying  assump¬ 
tion  that  A  =  0. 

ii.  By  inserting  the  softening  law  (3.13)  into  the  evolution  equation  (3.12)3  for  the 
softening  parameter  a,  noting  that  dq(j)  =  1  and  using  expression  (3.19),  one  obtains 

H~l  q  =  —a  =  —A  8s-  (3.22) 

Now,  the  stress-like  hardening  variable  q  must  remain  a  piecewise  smooth  function  if  the 
yield  criterion  <f>(r,q )  <  0  is  to  retain  its  classical  interpretation.  It  follows  that  the 
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preceding  relation  can  remain  mathematically  meaningful  only  if  the  softening  modulus 
itself  is  a  distribution,  i.e., 

H*1  =  H-1  6s  where  H  =  constant  <  0.  (3.23) 


This  result  is  amenable  to  a  compelling  physical  interpretation  which  relates  H  to  the 
fracture  energy  of  the  material,  see  Simo,  Oliver  &  Armero  [1993].  A  similar  interpretation 
holds  for  damage  models  (Oliver  &  Simo  [1994]). 

The  localization  condition  (3.20)  can  be  cast  into  the  following  equivalent  form 


LEMMA  3.2.  Let  Qe/  denote  the  perfectly  plastic  spatial  acoustic  tensor  evaluated  on  the 
surface  S  via  the  conventional  expression 


Q1  ■=  » 


c  — 


c dT<j>  ®  c dT<t> 


n  (for  A  >  0). 


dT<j> :  cdT(f>  J  5 

Then  the  localization  condition  (3.20)  is  equivalent  to  the  requirement 

QsPLvM  =  0  «=►  Ker[Q?]^0. 

which  holds  in  and  only  if  the  condition  detfQ^]  =  0  holds. 


(3.24) 


(3.25) 


PROOF.  Multiplying  both  sides  of  (3.20)  by  cdT<j)  and  contracting  the  result  with  dT<p 
yields  the  following  expression  for  the  multiplier:  A 

A  =  1  a -7  dr (t> :  c  [La [<p]  0  n\.  (3.26) 

dT<f>  :  cdr<p 

By  applying  the  tensor  c  to  both  sides  of  the  localization  condition  (3.20)  and  the  result 
to  the  normal  n  one  obtains 


[n  c  n  ]  La\}p\  —  A  [c  dT<f>]  n  =  0  (3.27) 

Inserting  (3.26)  into  (3.27)  and  rearranging  terms  yields  (3.25)  and  hence  the  result.  ■ 

It  is  worth  pointing  out  that  the  localization  condition  (3.25)  in  terms  of  the  acoustic 
tensor  does  not  involve  the  softening  parameter  H.  This  is  in  sharp  contrast  with  similar 
conditions  derived  in  the  literature  by  means  of  the  Thomas-Hill-Mandel  analysis  of  weak 
discontinuities  (acceleration  waves),  see  e.g.  Needleman  &  Tvergard  [1992].  An  identical 
condition  holds,  therefore,  for  the  perfectly  plastic  case  corresponding  to  Tt  =  0. 
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3.3  The  Limiting  Problem:  Plastic  Flow  on  the  Slip  Surface  S 

By  virtue  of  Assumption  3.1,  <£(S)  is  a  slip  surface  since  [y>]  •  n  =  0.  Our  goal  below 
is  to  identify  the  reduced  problem  that  governs  plastic  flow  on  the  slip  surface  S.  To  do 
so  observe  first  that  <5  is  a  material  surface  and,  therefore,  N  =  0.  As  a  result,  since 
n  -  F~tN ,  it  follows  that  L^n  =  0.  Time  differentiation  of  the  slip  condition  then 

yields  the  additional  relation 

0  =  *7-  ([y>]  •  n)  =  Ly  !v?J  •  n  +  M  •  L vn  =>  Ly  [y>J  •  n  =  0.  (3.28a) 

at 

Now  let  {ij,  f2,  n}  denote  the  convected  basis  at  y  =  <p(S),  so  that  ta  =  FY ,Q  is  a  normal 
parametrization  of  S.  Since  ‘f(S)  is  convected  by  the  flow  of  v  it  follows  that 

ta  •  n  =  0  for  a  =  l,2.  (3.286) 

The  localization  condition  (3.20)  together  with  (3.28a, b)  then  yield  the  following  results. 

i.  By  contracting  both  sides  of  (3.20)  with  n®n  and  ta  ®  tp,  a,  (3  =  1, 2  one  obtains 

n-[dT(f>]n  =  0  and  ta  ■  [dT<j>\tp  =  0  for  a,  ft  =  1,2,  (3.29a) 

while,  in  general, 

Ta  ’  [ dT<j)}n  7^  0  for  a  =  1,  2.  (3.296) 

Observe  that  for  J2-flow  theory  dT<p  =  dev[r],  where  dev[r]  stands  for  the  (Kirchhoff) 
stress  deviator.  Therefore,  for  J2-flow  theory  conditions  (3.29a, b)  assert  that  the  only  non¬ 
zero  stresses  on  a  slip— surface  are  the  resolved  shears  ra.  Motivated  by  this  interpretation, 
in  the  general  case  we  shall  refer  to  ra  defined  by  (3.29b)  as  the  generalized  resolved  shears. 

ii.  Since  S  is  material  and  ta  are  convected  by  it  follows  that 

4  (\<P\  ’  tot)  =  Lvl<p]  ■  ta  since  L^ta  =  0.  (3.30) 

at 

Therefore,  by  contracting  both  sides  of  (3.20)  with  n  ®  £a,  a  —  1,2  and  using  results 
(3.28a, b)  together  with  (3.30)  one  obtains 

■^AM-ta)  =  LvM-ta  =  2\Ta,  a  =  1,2.  (3.31) 

at 

It  is  clear  that  =  [<£>]  •  ta  represent  the  components  of  the  plastic  slip  taking  place  on 
the  surface  S.  Moreover,  since  4>(t)  is  homogeneous  of  degree  one,  the  yield  criterion  on 
S  can  be  written  as 


<f>(r, ?)|5  =  :=  \/Ti2  +  r22  _  +  ryl  -  0> 


(3.32) 
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As  a  result,  the  evolution  equations  (3.12)2,3,4  collapse  to  the  following  problem  posed  on 
the  slip  surface  S: 

£or  =  2A  Ta  1 

a  =  A  j  (3.33) 

A  >  0,  <t>s(ra,  a)  <  0  and  A  <f>s(ra, <*)  =  0. 

It  is  interesting  to  observe  that  this  problem  is  the  mathematical  transcription  of  the 
classical  Schmidt  law,  which  asserts  that  plastic  flow  on  a  slip  surface  is  proportional 
to  the  resolved  shear  stress  for  J2~ flow  theory.  See  Asaro  [1983]  for  a  micromechanical 
motivation  of  Schmidt  law. 

REMARK  3.4.  The  numerical  implementation  of  the  distributional  character  of  the  hard¬ 
ening  law  can  be  achieved  in  two  alternative  ways.  First  one  can  work  directly  with  the 
standard  problem  (3.12)  or,  equivalently,  with  the  model  of  multiplicative  plasticity  de¬ 
scribed  Remark  3.3.  Effective  numerical  algorithms  for  the  latter  form  of  the  model  are 
described  in  Simo  [1992].  The  distributional  character  of  the  softening  law  is  incorporated 
into  the  algorithm  by  replacing  the  delta  function  8s  in  the  softening  law  (3.13)  with  a 
delta-sequence,  as  described  in  the  next  section.  Alternatively,  one  can  work  directly  with 
the  limiting  problem  (3.33)  thus  by-passing  the  use  of  6-sequences.  Both  approaches  yield 
excellent  results  which  are  totally  mesh  independent.  The  former  approach,  however,  has 
the  advantage  of  involving  only  a  trivial  modification  of  existing  algorithms.  ■ 

4  A  New  Finite  Element  Method  for  Localization 

The  boundary  value  problem  to  be  solved  numerically  consists  of  the  weak  form  (2.11) 
of  the  equilibrium  equations,  with  test  functions  lying  in  the  space  V  defined  by  (2.12), 
supplemented  by  the  constitutive  model  described  above.  A  key  condition  to  be  satisfied 
by  the  spatial  discretization  emanates  from  the  analysis  and  conclusions  summarized  in 
the  preceding  section,  if  the  features  exhibited  by  the  continuum  solution  are  to  be  ac¬ 
curately  reproduced  by  the  numerics.  The  discretization  must  be  able  to  replicate  strain 
fields  that  exhibit  bounded  measures  (i.e.,  delta  functions)  in  order  to  capture  both  the 
localized  plastic  deformation  and  the  distributional  character  of  the  softening  law  that 
render  mesh-insensitive  numerical  solutions.  It  is  well-known  that  he  standard  Galerkin 
finite  element  method  cannot  meet  such  a  condition,  thus  leading  to  overly  diffused  re¬ 
sults  that  exhibit  a  strong  mesh  dependence.  The  technique  described  below  falls  within 
the  class  of  Assumed  Enhanced  Strain  (AES)  methods  proposed  in  Simo  k  Rifai  [1990] 
and  subsequenty  generalized  to  the  finite  deformation  regime  in  Simo  k  Armero  [1992] 
and  Simo,  Armero  &  Taylor  [1993].  Within  this  framework,  the  preceding  requirement  is 
achieved  via  a  particular  local  enrichment  of  the  deformation  field. 


A  New  Approach  to  Strain  Softening  in  Solids 


16 


(A)  Boundary  3£2h  =  (B)  Element  with  Regularization  Band 


FIGURE  4.1.  (A)  Finite  Element  discretization.  Definition  of  the  subdomain 
fth  with  boundary  dQh.  (B)  A  typical  element,  with  constant  normal  Ne 
to  the  element  shock  and  local  domain  involved  in  the  construction 
delt  a-sequences . 


4.1  A  Class  of  Assumed  Enhanced  Strain  Methods 

Consider  a  standard  finite  element  discretization  ft  &  U^Lj  Qe  of  the  domain  Q  with 
characteristic  mesh— size  h  >  0,  as  shown  in  Figure  4.1.  We  introduce  at  the  outset  the 
decomposition  (2.6)  of  the  displacement  field,  with  associated  test  functions  lying  in  the 
space  V  defined  by  (2.12).  For  convenience,  we  shall  denote  by  V  the  space  spanned  by 
the  regular  part  fj  of  the  test  functions  rj  —  fj  +  Ms/3  G  V . 

The  key  idea  in  the  method  described  below  is  to  work  directly  in  terms  of  the  regular 
part  of  the  displacement  field  and  the  admissible  test  functions,  which  are  amenable  to 
a  conventional  C° -finite  element  approximation  V^1  C  V.  For  the  sake  of  concreteness, 
consider  an  approximation  via  constant  strain  triangles  of  tetrahedra,  i.e, 

Vh  := 

The  singular  (distributional)  character  of  the  strain  field  in  the  continuum  problem  is 
captured  within  the  proposed  numerical  approach  via  a  local  enhancement  at  the  element 
level;  i.e.,  by  considering  enhanced  deformation  gradients  F ^  with  associated  variations 
denoted  by  Gh  of  the  form 


t  €  [c°m 


tt-dir 


fj e  G  (f2e)]ndim  and 


-h 

V 


=  0 


(4.1) 


F^  =  GRAD[y?ft]  +  j Feh  and  Ghe  =  GRAD  [77^]  -f-  G^  .  (4.2) 

Galerkin  Enhanced  Galerkin  Enhanced 

We  denote  by  £%  the  finite  element  spaces  of  enhanced  strains  Fh  and  let  £*  be  the  cor¬ 
responding  weighting  space  of  admissible  variations.  In  general  we  may  have 
a  situation  reminiscent  of  the  so-called  Petrov— Galerkin  formulations.  A  key  additional 
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condition  placed  on  both  and  £^  is  that  the  finite  element  approximations  be  discon¬ 
tinuous  between  elements.  Within  this  context,  the  appropriate  variational  formulation  of 
the  discretized  finite  element  problem  is  given  by  the  following  system: 


J  Ph  :  GRAD[»7ft]  dQ  =  J  f-i]hdfi  +  j  t  ■  T)h  <£T,  Vrjh  €  Vh,  j 

/  Ph:Ghedf2  =  0,  and  e  =  1,2,  •••,£?. 

Ji2e 


(4.3) 


The  nominal  stress  field  in  (4.3)  is  computed  via  the  specific  constitutive  model  under 
consideration,  but  evaluated  with  the  enhanced  strain  field  defined  by  (4.2)i.  The  design 
of  the  weighting  space  S!^  is  restricted  by  the  following  key  conditions  set  forth  in  Simo  & 
Armero  [1992],,  which  ensure  stability  and  consistency  of  the  of  the  AES  method. 

i.  Consistency:  For  the  linear  interpolations  (4.1),  strains  in  £^  must  have  zero  element 
mean,  i.e.,  Jn^  G\  dfi  =  0  for  all  Gh  €  £^. 

ii.  Stability :  The  space  GRAD[Vfc]  of  strains  associated  with  test  functions  in  Vk  and  the 
space  £if  must  have  null  intersection;  i.e,  GRAD[V/i]  fl  £i(  —  0. 


For  a  convergence  proof  restricted  to  linear  problems  see  Reddy  h  Simo  [1992]. 


4.2  Design  of  the  Enhanced  Strain  Fields 

To  describe  the  proposed  construction,  it  is  convenient  to  suppose  for  the  moment 
that  the  material  surface  S  is  given.  Furthermore,  for  simplicity  we  restrict  our  subsequent 
developments  to  the  case  ;m  =  2.  The  case  ndim  =  1  is  described  in  detail  in  Simo,  Oliver 
&  Armero  [1993].  Consistent  with  the  linear  interpolation  in  (4.1),  we  assume  that  the 
approximation  S ^  to  S  is  a  polygonal  line,  with  vertices  lying  on  the  sides  of  the  triangular 
elements,  which  intersect  an  element  Qt  along  the  segment  with  constant  unit  normal 
Ne.  Let  J  denote  the  set  of  numbers  associated  with  elements  that  intersect  Sh,  i.e., 

J  :=  (e€  {1,2,  ...,£}:  £2ef)Sh^&}.  (4.4) 

Now  let  N  be  the  global  node  numbers  associated  with  the  elements  Qe  intersecting  Sh; 
i.e.,  such  that  e  €  J.  We  set  N  =  N+  U  iV_,  where  N+  (respectively  iV_)  are  the  global 
node  numbers  in  N  associated  with  nodes  lying  ahead  of  (behind  of)  the  shock  S  .  Clearly, 
the  polygonal  lines  with  nodal  points  associated  with  the  sets  N±  define  approximations 
(also  denoted  by  <Si)  to  the  lines  Si  introduced  in  Section  2,  see  Figure  2.  With  this 
notation  at  hand  we  proceed  as  follows. 

(A)  Construction  of  the  space  €£.  Assuming  that  the  jump  [9?]  in  the  motion  is 
constant  within  a  typical  element  f?e,  e  €  J ,  and  in  view  of  the  decomposition  (2.6)  for 
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FIGURE  4.2.  Restriction  of  the  function  Ms  to  a  typical  element  Qe 

the  displacement  field,  it  is  natural  to  adopt  the  following  approximation  for  the  enhanced 
strain  field: 

£*  :=  | Fh  6  [L2{Q)\2  x  [L2(ft)]2  :  =  ae  ®  GRAD [M^J  with  a£e«2}.  (4.5) 

Here  ote  are  element  parameters,  which  are  suppose  to  replicate  the  jump  within  the 
element.  In  addition,  AfJ  is  the  restriction  to  the  element  Qe  of  an  approximation  to 
the  function  Ms  that  satisfies  the  two  conditions  i  and  ii  set  forth  in  Section  2.2  and  is 
constructed  as  follows. 

For  a  typical  element  Qe  with  e  G  J,Sh  intersect  two  sides  of  Qc  with  common  node 
which  will  be  denoted  in  what  follows  by  X*.  Let  N*  denote  the  element  shape  function 
associated  with  this  node.  Denoting  by  he  >  0  the  distance  from  X*  to  the  opposite  side, 
with  unit  outward  normal  Me  (see  Figure  4.1),  we  have 

N*{X)  :=  1  —  (X  —  X*)  • Me/he  with  \Me\  =  1.  (4.6) 

The  function  M<~  is  then  obtained  as  the  assembly  of  local  element  functions  MSe  defined 
as  (see  Figure  4.1) 

MSe  :=  HSe  -  N*e  where  F5e(X):=|j  x  G  S+’  (4'7) 

Here  f2e+  and  l?e_  are  the  regions  within  the  element  ahead  and  behind  of  the  shock.  This 
definition  clearly  ensures  that  |MsJ  =  1  (condition  i)  while  vanishing  outside  of  the  region 
region  Qh  limited  by  the  polygonals  and  the  boundary  of  the  body  (condition  ii).  The 
gradient  of  the  function  M$e  is  given  by  the  expression 

GRAD[MsJ  :=  8SeNe  -  Me, 


(4.8) 
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FIGURE  4.3.  Illustration  of  the  gradient  GRAD[Ms]  for  a  typical  element  fie 


which  is  illustrated  in  Figure  4.3. 


REMARK  4.1.  In  the  actual  implementation  of  (4.7)  the  H-function  Hse  is  replaced  by  an 
H-sequence.  Such  a  sequence  is  constructed  by  replacing  the  straight  segment  S'f  with  a 
narrow  band  of  thickness  k  denoted  by  12*.  The  derivative  of  the  function  MSe  regularized 
in  such  a  way,  which  is  denoted  by  M|e,  is  given  by  (see  Figure  4.4) 


GRAD[Mjj 


where  6$c 


i  jk  ifxei?*, 

0  otherwise  . 


(4.9) 


Observe  that  Me  =  Ne  only  if  the  shock  is  aligned  with  the  mesh.  The  regularization 
parameter  k  >  0  is  unrelated  to  the  characteristic  mesh  size  parameter  h  and  is  chosen  as 
small  as  possible  for  fixed ,  finite  h,  within  the  limitations  imposed  by  machine  precision. 
Within  this  context,  enforcement  of  the  distributional  character  of  the  hardening  law 
becomes  trivial:  The  delta-sequence  (4.9)2  induces  the  corresponding  delta-sequence  for 
the  (inverse  of  the)  softening  modulus 


kU  if 

oo  otherwise  . 


(4.10) 


The  value  7i^(X)  =  oo  for  X  $  Q*  is  consistent  with  elastic  behavior  outside  the  local¬ 
ization  band.  This  approach  is  totally  unrelated  to  the  so-called  characteristic  length,  an 
ad-hoc  concept  widely  used  in  smeared  crack  models,  see  e.g.,  Hilerborg  [1985],  Bazant 
[1986]  and  Oliver  [1989]  for  the  best  available  version.  ■ 

The  proposed  assumed  enhanced  strain  approach  is  clearly  strain  driven  and  fits, 
therefore,  the  usual  framework  of  constitutive  integration  algorithms  (see  Simo  [1992]  for 
a  comprehensive  review).  In  particular,  for  the  models  described  in  the  preceding  Section, 
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X,e£2* 


FIGURE  4.4.  Illustration  of  the  regularized  gradient  GRAD[MjJ  within  a 
typical  element  f2e- 


the  standard  return  mapping  algorithms  remain  unchanged;  the  only  modification  needed 
is  the  replacement  of  the  softening  modulus  with  the  sequence  (4.10). 

(B)  Construction  of  the  weighting  space  £k.  The  choice  adopted  below  is  motivated 
by  the  proof  of  part.(ii)  in  Lemma  2.1,  where  condition  (2.13b)2  is  obtained  essentially  by 
considering  in  the  weak  form  (2.11)  the  continuum  counterpart  of  the  enhanced  variations 
Gh  appearing  in  (4.2)2.  Relative  to  a  normal  parametrization  Xe  =  Ye  +  r]Ne,  where  Ye 
is  on  the  segment  and  —k  <  rj  <  k,  we  set 

£k  :=  {GA  €  [L\Q)f  x  [L\Q)f  :  Ghe  =  [SkSe  -  ^k']  {fie  ®  Ne) ,  ae  €  R2}  -  (4.11) 

The  function  ^{rj)  arbitrary  although  restricted  by  the  condition  ipk{r])U=k-^k(Tl)\v=-k  = 
1  which  ensures  satisfaction  of  the  stability  condition  on  -  Observe  carefully  that 
£*  ^  £b  except  in  the  case  Me  =  Ne  where  the  shock  is  aligned  with  one  triangle  side. 
The  ultimate  justification  for  (4.11)  lies  on  the  following  result  which  furnishes  the  finite 
element  counterpart  of  relation  (2.13)2: 


G\  :  Ph  dQ  =  0  VGh  <E  £k  4=4>  PkNe  - 


(4.12) 


Here  Pj  is  the  stress  on  the  band  Qke  and  Pf  is  the  stress  on  the  region  [Qe\£2k]  0  Qe+ 
ahead  of  the  shock.  The  proof  of  this  relation  is  straightforward  and  follows  from  the  fact 
that  Ph  is  piecewise  constant.  Observe  that  the  jump  condition  [PefeJiVe  =  0  is  proved 
exactly  as  in  Lemma  2.1.  Finally,  from  definition  (4.11)  it  is  immediate  to  verify  the 
£k  ft  GRAD[Vfe]  =  0.  Hence,  £k  defined  by  (4.11)  satisfies  the  consistency  and  the  stability 
requirements  for  AES  methods. 
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4.3.  A  Representative  Numerical  Simulation. 

To  illustrate  the  preceding  methodology  we  consider  the  isotropic  damage  model  de¬ 
scribed  in  Oliver  &  Simo  [1994]  in  an  uni-axial  tension  test.  It  can  be  shown  (see  the 
preceding  reference)  that  the  localization  condition  determines  the  angle  d  between  the 
normal  N  and  the  loading  direction  e  as  d  =  y/v,  where  v  is  the  Poisson’s  ratio.  A  typical 
numerical  result  for  v  =  0.4  is  shown  in  Figure  4.5  for  a  completely  non- structured  finite 
element  mesh.  The  computed  deformed  shape  of  the  specimen,  shown  in  Figure  4.5,  shows 
that  the  entire  deformation  takes  place  in  the  band  of  elements  intersecting  a  straight  line 
at  7r/2  +  d,  as  expected.  Outside  the  band  the  deformation  is  elastic  and,  therefore,  the 
specimen  behaves  nearly  as  two  rigid  solids  sliding  with  respect  to  each  other.  The  load 
displacement  curve,  not  shown  in  the  figure,  exhibits  absolutely  no  mesh  dependence. 


FIGURE  4.5.  Uniaxial  tension  test  for  am  isotropic  damage  model.  (A)  Non- 
structured  finite  element  mesh  and  computed  domain  12^  containing  the  shock 
Sh.  (B)  Computed  deformed  shape  exhibiting  the  strong  localization. 
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BING  C  YEN  /  IRVINE,  CA 

BRITISH  EMBASSY  /  ELLIS,  WASHINGTON,  DC 

BYU  /  ROLLINS,  PROVO,  UT 

CALIF  INST  OF  TECH  /  PASADENA,  CA 

CALTRANS  OFFICE  OF  RESEARCH  /  HOLLAND,  SACRAMENTO,  CA 
CATHOLIC  UNIV  /  CE  DEPT  (KIM)  WASHINGTON,  DC 
CENTRIC  ENGINEERING  SYSTEMS  INC  /  TAYLOR,  PALO  ALTO,  CA 
CHALMERS  UNIVERSITY  OF  TECHNOLOGY  /  TEPFERS,  412  74  GOTEBORG, 
CHEUNG  AND  ASSOCIATES  /  CHEUNG,  IRVINE,  CA 
COLORADO  SCHOOL  OF  MINES  /  GOLDEN,  CO 
COLORADO  ST  UNIV  /  FORT  COLLINS,  CO 
COMPUTATIONAL  MECHANICS  /  BREBBIA  SOUTHAMPTON, 

CORNELL  UNIV  /  ITHACA,  NY 

COUNTY  OF  VENTURA  /  TAKAHASHI,  VENTURA,  CA 
CRREL  /  KOVACS,  HANOVER,  NH 
CSU  CHICO  /  ARTHUR,  CHICO,  CA 
CSU  FULLERTON  /  RAMSAMOOJ,  FULLERTON,  CA 
DAMES  &  MOORE  /  LOS  ANGELES,  CA 

DET  NORSKE  VERITAS  RESEARCH  AS  /  BERGAN,  VERIT AS VEIEN  1,  N-1322  HOVIK 

DTIC  /  ALEXANDRIA,  VA 

ENVIROPLEX  /  AUDIBERT,  HOUSTON,  TX 

FAU  /  BOCA  RATON,  FL 

FAU  /  REDDY,  BOCA  RATON,  FL 

GEORGE  WASHINGTON  UNIV  /  ENGRG  &  APP  SCI  SCHL  (FOX),  WASHINGTON,  DC 
GEOTECHNICAL  ENGRS,  INC  /  MURDOCK,  WINCHESTER,  MA 
GEOTEST  INSTRUMENT  CORP  /  B ABEND REIER,  SILVER  SPRINGS,  MD 
GERWICK  INC  /  SAN  FRANCISCO,  CA 

GIANNOTTI  &  ASSOCIATES  OF  TEXAS  INC  /  GIANNOTTI,  VENTURA,  CA 

HEUZE  /  ALAMO,  CA 

HKS  INC  /  PAWTUCKET,  RI 

HQ  AFESC  /  TYNDALL  AFB,  FL 

IOWA  STATE  UNIV  /  AMES,  IA 

JOHN  HOPKINS  UNTV  /  COX,  BALTIMORE,  MD 

JOHN  HOPKINS  UNIV  /  LADE,  BALTIMORE,  MD 

KARAGOZIAN  AND  CASE  /  CRAWFORD,  GLENDALE,  CA 

LANTNAVFACENGCOM  /  CODE  411,  NORFOLK,  VA 


LAWRENCE  LIVERMORE  NATIONAL  LAB  /  MAKER,  LIVERMORE,  CA 

LAWRENCE  LIVERMORE  NATIONAL  LAB  /  MCCALLEN,  LIVERMORE,  CA 

LAWRENCE  TECH  UNIV  /  SOUTHFIELD,  ME 

LEHIGH  UNIV  /  BETHLEHAM,  PA 

LIN  OFFSHORE  ENGRG  /  SAN  FRANCISCO,  CA 

LOCKHEED  /  RSCH  LAB  (M.  JACOBY),  PALO  ALTO,  CA 

LOCKHEED  /  RSCH  LAB  (P  UNDERWOOD),  PALO  ALTO,  CA 

MAINE  MARITIME  ACADEMY  /  LIB,  CASTTNE,  ME 

MARC  ANALYSIS  RSCH  CORP  /  HSU,  PALO  ALTO,  CA 

MCCLELLAND  ENGRS  /  HOUSTON,  TX 

MICHIGAN  UNIV  /  HOUGHTON,  MI 

MIT  /  R.V.  WHITMAN,  CAMBRIDGE,  MA 

MOBILE  R&D  CORP  /  DALLAS,  TX 

MONTANA  STATE  UNIV  /  PERKINS,  BOZEMAN,  MT 

NATL  ACADEMY  OF  SCIENCES  /  NRC,  DR.  CHUNG,  WASHINGTON,  DC 

NAVFACENGCOM  /  CODE  04A3,  ALEXANDRIA,  VA 

NAVFACENGCOM  /  CODE  1002B,  ALEXANDRIA,  VA 

NAVFACENGCOM  /  CODE  163,  ALEXANDRIA,  VA 

NAVSHIP  /  CODE  245L,  FPO  AP, 

NCSC  /  PANAMA  CITY,  FL 

NEW  ENGLAND  MARINE  RESEARCH  LAB  /  LIB,  DUXBURY,  MA 
NFESC  ECDET  /  ESC61  (WU),  WASHINGTON,  DC 
NFESC  ECDET  /  ESC6I  CECILIO,  WASHINGTON,  DC 
NORDA  /  BENNETT,  NSTL,  MS 

NORTH  CAROLINA  STATE  UNIV  /  RAHMAN,  RALEIGH,  NC 
NORTHROP  /  CHEN,  HAWTHORNE,  CA 
NORTHWESTERN  UNIV  /  LIU,  EVANSTON,  IL 
NORTHWESTERN  UNIVERSITY  /  EVANSTON,  IL 
NORTHWESTERN  UNIVERSITY  /  BAZANT,  EVANSTON,  IL 
NORTHWESTERN  UNIVERSITY  /  CE  DEPT  (BELYTSCHKO),  EVANSTON,  IL 
NRL  /  VALENT,  STENNIS  SPACE  CENTER,  MS 
NSF  /  ASTILL,  ARLINGTON,  VA 

NSF  /  STRUC  &  BLDG  SYSTEMS  (KP  CHONG),  WASHINGTON,  DC 
NTH  /  LANGSETH,  N-7034  TRONDHEIM, 

NTH  /  MALO,  N-7034  TRONDHEIM, 

NUSC  /  LIB,  NEW  LONDON,  CT 

NY  MARITIME  COLLEGE  /  BRONX,  NY 

OCNR  /  CODE  10P4  (KOSTOFF),  ARLINGTON,  VA 

OCNR  /  CODE  1121  SILVA,  ARLINGTON,  VA 

ONR  /  CODE  1132SM,  ARLINGTON,  VA 

ONR  /  RAMBERG,  ARLINGTON,  VA 

OREGON  STATE  UNIV  /  CORVALLIS,  OR 

OREGON  STATE  UNIV  /  CORVALLIS,  OR 

OREGON  STATE  UNIV  /  CE  DEPT  (YIM),  CORVALLIS,  OR 

OREGON  STATE  UNIV  /  DEPT  OF  MECH  ENGRG  (SMITH),  CORVALLIS,  OR 

PENN  STATE  UNIV  /  LAB,  STATE  COLLEGE,  PA 

PMB  ENGRG  /  LUNDBERG,  SAN  FRANCISCO,  CA 

PORTLAND  STATE  UNIV  /  MIGLIORI,  PORTLAND,  OR 

PURDUE  UNIVERSITY  /  WEST  LAFAYETTE,  IN 

SAN  DIEGO  STATE  UNIV  /  SAN  DIEGO,  CA 

SCHWER  ENGR  &  CONSULTING  SERVICES  /  SCHWER,  BURR  RIDGE,  IL 
SCOPUS  TECHNOLOGY  INC  /  (B  NOUR-OMID),  EMERYVILLE,  CA 
SCOPUS  TECHNOLOGY  INC  /  (S.  NOUR-OMID),  EMERYVILLE,  CA 
SEAL  TEAM  6  /  NORFOLK,  VA 


SEATTLE  UNTV  /  SEATTLE,  WA 
SHANNON  AND  WILSON  INC  /  SEATTLE,  WA 
SJSU  /  VUKAZICH,  SAN  JOSE,  CA 

SOUTH  DAKOTA  SCHOOL  OF  MINES  AND  TECH  /  BANG,  RAPID  CITY,  SD 
SRI  INTL  /  ENGRG  MECH  DEPT  (GRANT),  MENLO  PARK,  CA 
SRI  INTL  /  ENGRG  MECH  DEPT  (SIMONS),  MENLO  PARK,  CA 
STANFORD  UNTV  /  APP  MECH  DIV  (HUGHES),  STANFORD,  CA 
STANFORD  UNTV  /  CE  DEPT  (PENSKY),  STANFORD,  CA 
STANFORD  UNTV  /  LAW,  STANFORD,  CA 
STATE  OF  CALIF  /  SACRAMENTO,  CA 

STRUCTURAL  ANALYSIS  PROGRAMS  INC  /  WILSON,  EL  CERRITO,  CA 

TEXAS  A&M  UNIV  /  COLLEGE  STATION,  TX 

TEXAS  A&M  UNIV  /  HERBICH,  COLLEGE  STATION,  TX 

TEXAS  A&M  UNTV  /  ROSCHKE,  COLLEGE  STATION,  TX 

THE  EARTH  TECH  CORP  /  ARULMOLI,  IRVINE,  CA 

TU  DELFT  /  DE  BORST,  2600  GA  DELFT, 

TU  DELFT  /  VAN  MIER,  2600  GA  DELFT, 

TUFTS  UNTV  /  SANAYEI,  MEDFORD,  MA 
UCLA  /  JU,  LOS  ANGELES,  CA 

UCSD  /  SCRIPPS  INST  OF  OCEANOGRAPHY,  LA  JOLLA,  CA 
UNTV  OF  CAL  BERKELEY  /  FILIPPOU,  BERKELEY,  CA 
UNIV  OF  CAL  BERKELEY  /  GOVINDJEE,  BERKELEY,  CA 
UNIV  OF  CALIF  BERKELEY  /  BERKELEY,  CA 

UNIV  OF  CALIFORNIA  /  MECH  ENGRG  DEPT  (BAYO),  SANTA  BARBARA,  CA 

UNIV  OF  CALIFORNIA  /  MECH  ENGRG  DEPT  (BRUCH),  SANTA  BARBARA,  CA 

UNTV  OF  CALIFORNIA  /  MECH  ENGRG  DEPT  (LECKIE),  SANTA  BARBARA,  CA 

UNIV  OF  CALIFORNIA  /  MECH  ENGRG  DEPT  (MCMEEKING),  SANTA  BARBARA,  CA 

UNIV  OF  CALIFORNIA  /  MECH  ENGRG  DEPT  (TULIN),  SANTA  BARBARA,  CA 

UNIV  OF  CALIFORNIA  DAVIS  /  CE  DEPT  (HERRMANN),  DAVIS,  CA 

UNIV  OF  CALIFORNIA  DAVIS  /  CE  DEPT  (KUTTER),  DAVIS,  CA 

UNTV  OF  CALIFORNIA  DAVIS  /  CE  DEPT  (RAMEY),  DAVIS,  CA 

UNTV  OF  CALIFORNIA  DAVIS  /  CTR  FOR  GEOTECH  MODEL  (IDRISS),  DAVIS,  CA 

UNIV  OF  COLORADO  /  CE  DEPT  (HON-YIM  KO),  BOULDER,  CO 

UNIV  OF  COLORADO  /  MECH  ENGRG  DEPT  (PARK),  BOULDER,  CO 

UNIV  OF  CONN  /  LEONARD,  STORRS,  CT 

UNTV  OF  CONN  /  LIBRARY,  GROTON,  CT 

UNTV  OF  DELAWARE  /  NEWARK,  DC 

UNIV  OF  HAWAII  /  HONOLULU,  HI 

UNIV  OF  HAWAII  /  KANEOHE  BAY,  HI 

UNIV  OF  HAWAII  /  HONOLULU,  HI 

UNIV  OF  HAWAII  /  HONOLULU,  HI 

UNTV  OF  ILLINOIS  /  CE  LAB  (PECKNOLD),  URBANA,  IL 

UNIV  OF  MICH  /  ANN  ARBOR,  ME 

UNIV  OF  N  CAROLINA  /  CE  DEPT  (GUPTA),  RALEIGH,  NC 

UNIV  OF  N  CAROLINA  /  CE  DEPT  (TUNG),  RALEIGH,  NC 

UNIV  OF  NY  /  BUFFALO,  NY 

UNIV  OF  RHODE  ISLAND  /  KOVACS,  KINGSTON,  RI 

UNIV  OF  RHODE  ISLAND  /  VEYERA,  KINGSTON,  RI 

UNTV  OF  TENNESSEE  /  KNOXVILLE,  TN 

UNIV  OF  TEXAS  /  JIRSA,  AUSTIN,  TX 

UNTV  OF  TEXAS  /  TASSOULAS,  AUSTIN,  TX 

UNIV  OF  WASHINGTON  /  MATTOCK,  SEATTLE,  WA 

UNIV  OF  WISCONSIN  /  MILWAUKEE,  WI 

UNTV  OF  WYOMING  /  CIVIL  ENGRG  DEPT,  LARAMIE,  WY 


USA  BELVOIR  /  FORT  BELVOIR,  VA 
USACOE  /  BRADLEY,  VICKSBURG,  MS 

VA  POLY  INST  AND  STATE  UNTV  /  MITCHELL,  BLACKSBURG,  VA 
VIRGINIA  INST  OF  MARINE  SCI  /  GLOUCESTER  POINT,  VA 
WEIDLINGER  ASSOCIATES  /  LEVINE,  LOS  ALTOS,  CA 
WEST  VIRGINIA  UNIV  /  BARBERO,  MORGANTOWN,  WV 
WEST  VIRGINIA  UNIV  /  KIGER,  MORGANTOWN,  WV 
WOODS  HOLE  OCEANOGRAPHIC  /  LIB,  WOODS  HOLE,  MA 
WOODWARD  CLYDE  CONSULTANTS  /  OAKLAND,  CA 
WORCHESTER  POLYTECH  /  SULLIVAN,  WORCESTER,  MA 


